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This paper presents a framework for calibrating computational mod- 
els using data from several and possibly dissimilar validation experiments. 

The offset between model predictions and observations, which might be 
caused by measurement noise, model-form uncertainty, and numerical er- 
ror, drives the process by which uncertainty in the models parameters is 
characterized. The resulting description of uncertainty along with the com- 
putational model constitute a predictor model. Two types of predictor mod- 
els are studied: Interval Predictor Models (IPMs) and Random Predictor 
Models (RPMs). IPMs use sets to characterize uncertainty, whereas RPMs 
use random vectors. The propagation of a set through a model makes the 
response an interval valued function of the state, whereas the propagation 
of a random vector yields a random process. Optimization-based strategies 
for calculating both types of predictor models are proposed. Whereas the 
formulations used to calculate IPMs target solutions leading to the inter- 
val value function of minimal spread containing all observations, those for 
RPMs seek to maximize the models’ ability to reproduce the distribution 
of observations. Regarding RPMs, we choose a structure for the random 
vector (i.e., the assignment of probability to points in the parameter space) 
solely dependent on the prediction error. As such, the probabilistic de- 
scription of uncertainty is not a subjective assignment of belief, nor is it 
expected to asymptotically converge to a fixed value, but instead it casts 
the model’s ability to reproduce the experimental data. This framework 
enables evaluating the spread and distribution of the predicted response of 
target applications depending on the same parameters beyond the valida- 
tion domain. 
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Acronyms 


CDF Cumulative Distribution Function 

DGM Data Generating Mechanism 

IID Independent and Identically Distributed 

IPM Interval Predictor Model 

LS Least Squares 

PDF Probability Density Function 

RPM Random Predictor Model 


I. Introduction 

Model calibration refers to the process of prescribing the values of the parameters of 
a computational model according to experimental observations. Model-form uncertainty, 
measurement noise, and numerical error often prevent us from conhdently prescribing a 
hxed constant value for such parameters. Consequently, a family of parameter values is 
prescribed such that the predictions resulting from evaluating the computational model at 
any member of this family are a sufficiently accurate representation of the observations [4]. 

Several model calibration techniques are available in the literature. Many of them assume 
the structure 

y = M{x,p)+T], (1) 

where ?/ G M is the response or model’s output, M is the computational model, x G M"" 
is the state or explanatory variable, p G is the model’s parameter, and r; G M is a 
random variation caused by noise and measurement error. Many of the model calibration 
techniques are based on this model structure, the assumption of p being a hxed but unknown 
constant (i.e., the uncertainty in p is epistemic), and the assumption of p being independent 
and identically distributed (IID) following a normal distribution with zero mean and a hxed 
variance. In contrast, the framework proposed herein is applicable to models M, where p 
can either be epistemic or aleatory (i.e., uncertainty caused by inherent, and irreducible 
variability) taking on an arbitrary form*^. 

A typical regression problem consists of estimating the value of p in M given the set of ob- 
servations (xj, Pi), i = 1, . . . , A^, where N > Up. This is often carried out by searching for the 
parameter realization that minimizes the sum of squared residuals, pls- The precision of this 
estimate, which prescribes how much it can deviate from its “true value” within an epistemic 
framework, is often evaluated using conhdence intervals [8, 5]. The calculation of conhdence 
intervals, prediction intervals (i.e., intervals where future observations are expected to fall) 
and credible intervals [4] require having a probabilistic description of the uncertainty in p. 
This description often requires (i) knowing the distribution for the prediction errors, (ii) M 
and p taking particular forms (e.g., M depends linearly on p and the noise p is additive) or 
(iii) M being accurately represented by a linear approximation. As such, the suitability of 
the resulting conhdence and prediction intervals depends tigthly on the validity and accuracy 
of the underlying assumptions and approximations. 

^Equation (1) can be seen as a particular case, in which all epistemic uncertainties enter into M, whereas 
r] is a single additive aleatory uncertainty having a known structure. 
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The most common approach to model calibration is Bayesian inference. In parametric 
Bayesian inference the objective is to describe the model’s parameters as a vector of possibly 
dependent random variables using Bayes’ rule [4], The resulting vector, called the posterior, 
depends on an assumed prior random vector, and the likelihood function; which in turn 
depends on the set of N observations available and on the model structure of M. Whereas 
this approach does not make any limiting assumptions on the manner in which M depends 
on p, nor on the structure of the resulting posterior; it requires that p be epistemic. The 
presence of aleatory uncertainty and model-form uncertainty yields uncertainty character- 
izations for p that fail to describe the prediction error. This can be mitigated by adding 
a non-deterministic discrepancy term to M [4]. This practice however, prevents rolling-up 
the resulting prescription of uncertainty through other computational models. By roll-up we 
mean the process of calibrating the parameters of one or several computational models and 
using the resulting characterization of p to make predictions of a target application which 
also depends on p. In spite of this limitation, its high computational demands, and of the 
potentially high sensitivity of the posterior to the assumed prior, this method is commonly 
regarded as the benchmark for model calibration. 

This paper presents techniques for characterizing the uncertainty in parameters common 
to several computational models based on multiple validation experiments. These experi- 
ments might measure dissimilar quantities of interest having different state variables ranging 
over different domains (e.g., some might be temperatures at particular locations whereas 
others might be heat flux at other locations). For each of these validation experiments we 
have a computational model, and a set of observations. Uncertainty is characterized accord- 
ing to the dispersion of the prediction errors. This characterization and a computational 
model constitute a predictor model. This paper focuses on Interval Predictor Models (IPM), 
for which uncertainty is characterized as a hyper-rectangular set’’, and Random Predictor 
Models (RPM), for which uncertainty is characterized as a random vector supported in a 
hyper-rectangular set. Hyper-rectangular sets are preferred since they have the advantage 
that the range of an individual parameter is not affected by the value taken by the others. 
IPMs and RPMs are used to generate informative predictions of the spread of the response at 
values of the state that might extend beyond the range of observations (i.e., extrapolation). 
Such predictions are solely based on the structure of the computational models and their 
ability to reproduce the observations. 

The calibration of models that depend linearly in p and polynomial in x can be done 
rigorously [3]. This process yields a formal description of the model’s reliability, i.e., the 
probability that future observations will lie within the predicted range of responses by an- 
alytical, non-statistical means; of the uncertainty in p, and of the spread in the predicted 
response y. This paper extends these ideas to cases in which (i) data from multiple validation 
experiments are available (i.e., > 1, where is the number of validation experiments), 

(ii) the model structure depends arbitrarily on parameters and states and might not even 
have an explicit form (e.g., M is a hnite element model predicting the stress at a particular 
location of a nonlinear structure), and (iii) the resulting characterization of uncertainty is 
probabilistic. This generality precludes the formal descriptions listed above^. 

'^By hyper-rectangular set, we mean a set which can be represented by a Cartesian product of intervals. 
This constrains its orientation so that its edges are parallel to the coordinate axes. 

'^The inability to assess the model’s reliability precludes making any statement about how the predicted 
range of responses relates to future data. As such, the informative character of the prediction is limited. 
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This paper is organized as follows. Section II introduces basic notions and definitions and 
describes the main objectives of this article. This is followed by Section III where a framework 
for calculating IPMs is presented. Section IV introduces RPMs and a formulation for their 
calculation. This is followed by Section V where the methodology is extended to the roll-up 
problem. Finally, some conclusions are made. 

II. Problem Statement 

In the developments that follow we assume that experimental observations are stationary 
and independent, that the response of the models depend continuously on both the param- 
eters and the states, and that the models have a deterministic mathematical structure, i.e., 
they yield a fixed value of the response for fixed values of the state x and of the parameter p. 
Consider the ensemble^ of computational models = M^{x^ ,p), for j = 1, . . . n^, where the 
response is a scalar, the state x^ is a vector of length Uj, p ^ is a parameter, and Uy 
is the number of validation experiments. Note that all riy computational models share the 
same parameter p. The values of parameters that are not common to all Uy models can be 
set to fixed constant values, e.g., the corresponding term in the least squares solution, before 
applying this methodology. This practice enables rolling up the uncertainty in all calibrated 
models to system-level target applications. 

Further assume we have experimental observations corresponding to each of the Uy val- 
idation experiments. Denote by the number of observations corresponding to the jth 
validation experiment. The Data-Generating Mechanism (DGM) associated with each of 
these experiments is postulated to act on the state to produce a response. The responses 
can depend on state variables and on some other influences such as intrinsic variability and 
measurement noise. Let C M”'’ be a set of state variables, and C M be the set of 
responses corresponding to the DGM governing the jth validation experiment acting on ele- 
ments of VL The data sequence associated with this experiment is denoted as = (xl,yl), 
where i = 1, . . . iVL The DGMs might depend on parameters and states that are absent 
from the computational models. As compared to a DGM, the model might be subject to 
parametric and model-form uncertainty. Parametric uncertainty can be epistemic, aleatory 
or a combination of both. As compared to a computational model with the highest fidelity 
possible, often suffers from approximation, numerical, and discretization errors. These 
errors are a consequence of making computationally tractable. 

The first objective of this paper is to prescribe p according to the data sequences z^, z^, . . ., 
z””, based on the predictive ability of the computational models M^, . . ., M"” describing 

such experiments. The second objective is to use this prescription to generate informative 
predictions corresponding to unobserved realizations of the state. An informative prediction 
can be interpreted as a prediction that is consistent with salient features of the data sequences 
and with the predictive ability of the Uy computational models used to characterize p. In 
the context of IPMs, a prediction is an interval of possible responses whereas for RPMs a 
prediction is a random variable. The spread and probability concentration of such predictions 

‘^The following conventions will be used regarding the usage of subindices and superindices on the model’s 
state and model’s response. Subindices will be used to describe experimental realizations whereas su- 
perindices will be used to describe different states/responses; e.g., is the third observed response associated 
to model whose state and response are x'^ and respectively. 
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are driven by the model’s ability to reproduce the observations. As such, models that fail 
to properly capture the physics driving the DGM will yield large uncertainties and wide 
prediction ranges. 

For simplicity in the presentation, the two sections on IPMs and RPMs that follow assume 
that there is only one validation experiment (i.e., = 1). Consequently, we will drop the 

super-indices from the notation of states and responses. 

III. Interval Predictor Models 

An IPM is simply a function that returns an interval as output. These models were 
proposed in the framework of differential inclusions and set- valued dynamical systems [1, 
2], and further developed for the case in which prior knowledge on the data-generating 
mechanism is available [6, 7]. In the context of this paper, an IPM is a rule that assigns to 
each instance vector x & X a. corresponding outcome interval in Y . That is, an IPM is a 
set- valued map 

/ : X — Iy{x) C y, (2) 

where x is a state on which the model’s output depends, and Iy{x) is the prediction interval. 
Let M be any functional acting on a vector x of state variables and a vector p of parameters 
to produce an output y; i.e., y = M{x,p). A parametric IPM is obtained by associating to 
each X G A the set of all possible outputs y that result from varying p over some subset P 
of parameter space: 

Iy{x,P) = {y ■ y = M{x,p) for all p g P}. (3) 

Iy{x,P) will be an interval as long as M{x,p) is a continuous function of x and p, and 
P is a connected set. All instances of M and P considered in this paper will satisfy these 
restrictions. Attention will be limited to the case in which P is the bounded® hyper-rectangle: 
P = {p : pmin < P < Pmax}- The prediction of the IPM can be described as Iy{x,P) = 
[y{x,P), y(x,P)], where 

p(x, P) = min{M(x,p)}, (4) 

~ p£P 

p(x, P) = max{M(x,p)}. (5) 

p£P 

The functions y and y are envelopes of the family of inhnitely many predictions that result 
from evaluating M(x,p) for each parameter realization p & P. Note that the envelopes 
themselves might not be members of the family (e.g., the envelopes of a family of polynomials 
need not to be polynomials). The spread of ly, which is the separation between the envelopes, 

is 

6y{x,P) =y{x,P) -y{x,P). ( 6 ) 

The solutions to Equations (4-5) occur at vertices of P when M is a linear function of p 
(the particular vertex might vary with x). Otherwise, the solution might correspond to an 
interior point of P. The case in which M is a convex function of p can be efficiently solved 
by using convex optimization techniques (when solving for y) and by evaluating M at the 
vertices’s of P (when solving for y). If none of the above cases hold, nonlinear optimization 
or sampling techniques can be readily used. 

^Vector inequalities hold component-wise. 
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The DGM is often approximated by the Least Squares (LS) prediction y = M(x,pls )5 
where the parameter estimate pls is given by 

{Vi- (7) 

Whereas the LS prediction describes the overall trend of the data as described by a point 
estimate of p, Iy{x) describes its spread. A set P, which prescribe such a spread, might not 
even contain pls- 

Formulations for calculating two types of IPMs are presented next. These formulations 
parameterize the geometry of P differently. Denote by 9 the set of parameters required to 
fully prescribe P. The two parameterizations to be considered are (i) 9 = {c,m,a}, where 
c is the geometric center, m is the unit diagonal vector, and a > 0 G M is the length of the 
semi-diagonal of P; and (ii) 9 = {pmin,Pmax}, where pmin and pmax are the “lower left” and 
“upper right” corners of P. The dependence among these variables is given by 

Pmax T Pmin Pmax Pmin 

C — , TTl jj r , CX 

^ llPmax Pmin||2 

A. Hyper-rectangular Sets with Fixed Parameters 

Here we seek the IPM given by (3) where the uncertainty set P(6*) is parameterized by 
9 = {c,m,a} according to 


llPmax Pmin II 2 


( 8 ) 


f N 

Pls = argmin ^ 
P [ i=i 


P{9) = {p : c — am < p < c + am}. (9) 

The following optimization program yields an IPM whose uncertainty set (9) has hxed given 
values for c and m whereas a is optimized. 

The optimization programs in this paper utilize Ej.[-], the expected value operator with 
respect to x, in their objective functions. When x has a known distribution, the implied 
integral can be evaluated analytically or numerically depending on the integrand. Otherwise, 
the sample mean of the integrand using the data in z can be used to approximate it. 

Optimization Program 1. Consider the IPM defined by (3) and (9). Assume that the values 
c and m > 0 are set in advance and denote by Ex[-] the expectation operator with respect to 
the state variable x. The value of a required to fully prescribe this IPM is 

a = argmin [E^[6y{x, P)] : y{xi, P) < yi < y{xi, P), 1 < f < A} . (10) 

o>0 

Therefore we search for the set P with a hxed geometric center and unit diagonal vector 
that minimizes the expected interval spread such that all the observed responses are within 
the limits of the interval valued function ly. The corresponding uncertainty set, denoted as 
P, is given by (9) evaluated a.t 9 = {c, m, a}. Note that if Pi C P 2 , then 6y{x, Pi) C 6y{x, P 2 ), 
and Esc[Sy{x, Pi)] < Ex[Sy{x, P 2 )]. Consequently, the minimization in (10) yields uncertainty 
sets of minimal size in this sense. 

A few comments regarding the selection of the hxed parameters in P are in order. A 
natural candidate for c is pls- Even though this practice leads to a prescription of uncertainty 
in which the prediction corresponding to the center of the box optimally describes the overall 

6 of 24 


American Institute of Aeronautics and Astronautics 


trend of the data, the corresponding spread in the prediction might be snboptimal (i.e., the 
selection of another center point leads to a smaller expected spread). Regarding the nnit 
diagonal vector, the components of m shonld be made proportional to the expected range of 
nncertainty in the corresponding parameters. Note that m shonld not only acconnt for snch 
ranges bnt also for the different nnits in the corresponding p’s. 

The above formnlation does not take into acconnt nncertainty in experimental observa- 
tions. Error dne to instrnment calibration, data acquisition, and variation in experimental 
conditions yield random and bias errors. This uncertainty is often prescribed as error bars 
centered about the measurements, e.g., the uncertainty in {xi, yi) is cast as [xi ± dx, yi ± dy) 
where dx > 0 and dy > 0. The effects of experimental uncertainty can be taken into account 
by replacing the inequality constraints in (10) with 

y{x,P) — dy < yi < y{x,P) + dy, for all x G [xj — dx, Xt + dx], (11) 

The particular case in which dx = 0 yields p(xj, P) + dy > y^ > y{xt, P) — dy. The 
general form of (11) is hard to evaluate since the inequality applies to a continuum of 
state values. A natural approximation to (11) is given by y{xl,P) < yl < y{xl,P) for 
k = 1, . . . , 2”^+^, where {(x^,p^) : 1 < enumerates the vertices of the box 

[xi — dx, Xi + dx]x[yi — dy, yi + dy]. The formulations herein ignore experimental uncertainty, 
thus dx = 0 and dy = 0. Their extension to the general case can be carried out by taking 
these considerations into account. 

The formulation in (10) as well as some of the concepts to be introduced later in the paper 
can be described using the Zero-Prediction-Error Manifold. This is the manifold of parameter 
points for which the computational model M can exactly reproduce an observation point. 
More specihcally, the Zero-Prediction-Error Manifold for the data point {xi,yi) is dehned as 

Si = {p : yi- M{xi,p) = 0}. (12) 

Note that a feasible solution to (10) leads to a set P whose intersection with Si, S 2 , . . . , 
S]\f are all non-empty. When M depends linearly on p such manifolds are hyperplanes. 
When the response yi is outside the range of predictions associated with the model at Xj, 
Si is empty. When the observation is within such a range, there might be inhnitely many 
parameter realizations and, thus, inhnitely many hxed parameter point predictions M{xi,p) 
passing through yi 

The separation between c and the manifold Si in (12) can be evaluated via the Parametric 
Safety Margin (PSM), which is dehned as 


Pi = min{||c-p||“ : yi = M{xi,p)}, 

V 


(13) 


where the m-inhnity norm of a vector a is dehned by 


|a||“ = max 


|®fc| 

rrik 


The PSM Pi is the length of the semi-diagonal of the smallest hyper-rectangle oriented as 
P (i.e., same center and same diagonal orientation) for which one of the corresponding IPM 
envelopes passes through the ith observation point. Hence, pi evaluates the extent by which 


7 of 24 


American Institute of Aeronautics and Astronautics 


the presence of the data point {xi,yi) contribntes to the spread of P. The spread of the 
resnlting prediction and of the nncertainty set are driven by the observation(s) attaining 
the largest PSM valnes. Note however, that the particnlar observations in z attaining snch 
valnes vary with c and m. Eqnation (13) is eqnivalent to 

Pi = min {a : a > 0, yt = M{p, Xi), —am < c — p < am} , (14) 

p^a 

which, as opposed to (13), does not snffer from derivative discontinnities. Note that a = 
maxj=i„.Ar{pi}. PSMs are nsed in [3] to identify and eliminate ontliers from the data set and 
to calcnlate IPMs having tighter predictions. 

1. Fixed-parameters IPM Example 

In this example we will consider Validation Experiment 1. This experiment consists of 
= 50 observations of an ontpnt which depends on a single state variable {n} = 1) and 
two parameters {up = 2). This DGM depends nonlinearly on the state and parameters, 
and it is snbject to state-dependent measnrement noise (see Appendix). Fignre 1 shows the 
elements of the corresponding data seqnence, z^. 
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Figure 1. Observations, LS prediction, and confidence band. 

A mathematical model M^{x^,p) describing this DGM was bnilt. This model, which de- 
pends nonlinearly on G M and p G snffers from model-form and epistemic nncertainty. 
The prediction corresponding to the LS parameter, which is plsi = [0.9420, 1.0689]^, and 
the prediction interval fnnction [5] associated with the nonlinear regression are also shown 
in Fignre 1. In statistical inference [8], a prediction interval is an estimate of an interval 
in which fntnre observations will fall, with a certain probability, given what has already 
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been observed. Note that more than 60% of the observations lay outside the interval valued 
function. 

The model and the = 50 observations were used to build an IPM. The 

resulting predictor model will be called IPMl. The solution to (10) yields a = 0.2625, and 
Ex[5y] = 2.2505. Figure 2 shows the envelopes of IPMl for c = pls and m = [\/(h5, \/(h5]''^. 
All 50 observations are within the envelopes and the lower envelope passes through one of 
them, indicating that the corresponding constraint is binding. This particular observation 
point, located near (0.7,1.25), ends up prescribing the extension of P. Recall that the 
functional form of the envelopes depends exclusively on the observations and the structure 
of M^. Whereas M^(x^,pls) describes the dominant trend of the data by weighing all data 
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Figure 2. IPM with fixed parameters. 

points equally, the interval valued function Iyi{x^,P) describes their spread. This spread, 
which is driven by extreme observations, is not necessarily centered about M^(x^,pls)- In 
general, there is no basis to expect that (i) Pls ^ P, or that (ii) y^{x^) < M^(o;^,pls) < 
for all x^ G X, or that (hi) the centerline of Iyi{x^,P), {y^{x^,P) +y^{x^, P))/2, accurately 
describes the overall trend of the data. 

B. Hyper-rectangular Sets with Free Parameters 

In this formulation we seek the IPM in (3) for the uncertainty set P{0) parameterized by 

0 = {Pmin,Pmax} aCCOrdiug tO 


P(6I) = {p : < p < Pmax}, 

such that all components of 6 are optimal. 


( 15 ) 
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Optimization Program 2. Consider the IPM defined by (3) and (15). This IPM is fully 
prescribed by the limits of P, which are given by 

{pmm.,pmax) = ?vrgm.m{EjPy{x, Q)] : Q) < Vi < y{xi, Q), I <i < N, (16) 

Pa,Pb 

Q = p[0) where 0 = {pa,Pb}, Pa < Pb}- 

Therefore we search for the limits of P that minimize the expected interval spread such 
that all the observed responses are within the limits of the interval valued function Iy{x). 
The resulting uncertainty set, P, is given by (15) evaluated at 0 = {pmax, Pmin}. Remarks 
regarding the evaluation of E^[6] made in the previous section apply here as well. The PSM 
associated to resulting IPM are given by (13) with the parameters c and m corresponding 
to P. 

The additional flexibility of (16), as compared to (10), often yields IP Ms with tighter 
predictions and smaller uncertainty sets. As before, the spread of both the prediction and 
the uncertainty set are often driven by a few data points. These data points, which attain 
the largest PSM values from the full ensemble in z, often deviate signihcantly from the rest 
of the observations. The empirical Cumulative Distribution Function (CDF) of the PSMs 
is dehned to be the function whose value at an arbitrary real number p is the fraction of 
the N values of the PSMs which are less than or equal to p. Examining this empirical CDF 
enables visualizing how likely are the extreme observations prescribing the IPM envelopes. 
The majority of PSM values might be considerably smaller than the largest PSM value(s) 
of the set, even though the IPM is prescribed according to the data point (s) close to such 
values. 

1. Free-parameters IPM Example 

The free-parameter formulation in (16) was used to build an IPM based on the same compu- 
tational model and the same observations used in the previously. The resulting model, to be 
referred to as IPM2, is shown in Figure 3. The optimal solution leads to c = [0.8837, 0.9740]^, 
m = [0.3236,0.9462]^, and Ex[dy\ = 1.9642. The improved tightness in the model predic- 
tion of IPM2 as compared to that of IPMl is apparent. Note that most improvements are 
reflected in a tighter upper envelope, while the lower envelope degrades slightly. Further 
notice that the LS prediction is closer to upper envelope than in the IPMl case. 

Figure 4 shows the empirical CDF of the PSM associated with IPM2. The vertical lines 
are the 10-percentiles. Note that the PSM values range from 0.0005 to 0.3263, but about 
80% of them take on values that are less than 0.15. The observation(s) for which the PSM 
attains the largest value(s) prescribe the IPM whereas the rest of them are inconsequential. 
The optimization programs above can be cast as min-max problems since they lower the 
worst-case PSM (i.e., the 100-percentile of p) as much as possible. When outliers are present 
in the data set, or when the analyst wants tighter predictions for most of the observations 
(i.e., lower values of p), the minimization of lower percentiles can be considered instead [3]. 

IPMs are effective for describing extreme and possibly unlikely predictions but lack the 
hdelity to describe the likelihood of occurrence within lyilpjp). RPMs, to be introduced next, 
provide a means to mend for this dehciency. 
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Figure 3. IPM with free parameters. 



Figure 4. Empirical CDF of the PSMs associated to IPM2 
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IV. Random Predictor Models 


An RPM is simply a rule that assigns to each state vector x ^ X a. corresponding random 
variable in the output space Y. That is, an RPM is a random variable-valued map 

R : X ^ Ry{x) C Y, (17) 

where x is the state, and Ry{x) is a random process having a support that lies in Y. The 
prediction of the RPM is 

Ry{x) = {y - y = M{x,p), p ~ fp{p) for p e p(0)}, (18) 

where fp{p) is the joint Probability Density Function PDF of p having P{0) as its support set. 
As with the IP Ms presented above, we will focus on hyper-rectangular P’s parameterized by 
9. When both P and M are bounded, which is the case assumed in this paper, the support 
of Ry{x) is bounded as well. The PDF of the predicted response M{x,p) at any value of x 
is fully prescribed by that of p. The PDF of the random variable Ry{x) will be denoted as 
fM{x,p{ 0 )){y)- Therefore, statistics of the response, such as the mean Py{x) = Ep[M{x,p)], the 
variance i^y{x) = Ep[{M{x,p) — /ij^(x))^], and the support Iy{x) (given in Equations (4,5)); 
vary with x. Note that the limits of the support of an RPM are the envelopes of an IPM 
having the same P. 

Figure 5 shows an example of the RPM in (18). This hgure shows the LS prediction, the 
expected/mean response, and the 5-percentile curves for a uniform PDF supported in the P 
corresponding to IPM2. This RPM, which will be compared to other RPMs below, will be 
referred to as RPMl. Recall that IPM2 leads to a range of predictions with minimal spread. 
Note the sizable offset between the LS prediction and the expected response Py{x). It turns 
out that the mean response underestimates 41 of the 50 experimental observations by a 
considerable margin, thus, it is a poor descriptor of the overall trend of the data. Further 
notice that the probability of occurrence of the data sequence z under the IID assumption, 
which is the product of the individual probabilities at the observation points and whose 
values can be inferred from the percentile curves, is comparatively low (see Figure 8 for 
comparison) . 

This paper proposes a process for constructing a probabilistic description for p, fp{p) with 
p G P{9), based on a non-subjective hgure of merit. This hgure of merit rewards parameter 
realizations (i.e., assigns a higher value) leading to model predictions that closely describe 
the experimental observations (i.e., points in the parameter space lying on the manifolds S 
in (12)). The closer the parameter point to the manifolds associated to all the experimental 
observations, the larger the value of the hgure of merit, thus, the larger the value of the joint 
density function at such a point. Note that the “distance” from any parameter point to any of 
such manifolds, which can be evaluated using several diherent norms, will only depend on the 
experimental observation [xi,yi) and on the ability of the computational model y = M{x,p) 
to replicate such an observation. Note that these two underlying criteria also prescribe the 
likelihood function used in Bayesian inference. The developments that follow propose a few 
norms to evaluate such a “distance”, thus, a hgure of merit, and describe a procedure from 
which probabilistic descriptions of p, consistent with this hgure of merit, can be built. 

Assume the joint PDF for p is given by 

fp(p) = ^ foi’ P e (19) 

Go 
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Validation Experiment 1 



Figure 5. RPM for a uniform joint PDF (RPMl). 


where Cq is a normalization constant so the volume under the joint PDF over P{9) is one, and 
6i is a variable to be determined. The function 7(p) assigns a hgure of merit to each possible 
realization of p. Denote by d{p, S',) > 0 a norm of the separation between the parameter 
point p and the zero prediction error manifold in (12). The structure of the function 7 to be 
considered is given by 

, ( 20 ) 

where 0 < /?i 1 is a small number used to prevent 7 from being unbounded, and /32 > 0 

is a scaling parameter used to control the smoothness of 7(p). The closer p is to parameter 
points where the computational model is able to reproduce all experimental observations, 
the larger the value of 7, thus, the larger the value of at p. Note that in an idealized 
case in which there is no measurement noise or model-form uncertainty, and the only source 
of uncertainty is epistemic, 7(p) will be maximal at the “true” value(s) of p. 

Several norms d{p, Si) can be used to fully prescribe (20). Each of them will lead to a 
different joint PDF, thus to a different random process Ry{x). Three suitable norms are 


d{p,Si) = min{||p-p||2 : Pi = M{xi,p)}, 

p 

(21) 

d{p,Si,m) = pi = 

min{||p - p||~ : pi = M{xi,p)], 

p 

(22) 

d{p, a: 

{pi- M{p,Xi))‘^ 
Ci 

(23) 


N 


7(p) = A + /^2 5]] 


2=1 
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where Q is a normalization constant used to ensure that d{p, Si) < 1 for alH = 1 . . . iV. The 
norms in (21) and (22) measure the separation between any parameter point p and the zero- 
prediction-error manifold by evaluating a distance in the parameter space. Conversely, (23) 
measures separation by evaluating a distance in the output space. As such, the norm in (23) 
yields results that are structurally different from those based on the other two norms (unless 
M is linear in p). Furthermore, the norms in (21) and (22) require solving an optimization 
problem, whereas (23) only requires a single model evaluation. This makes the norm in 
(23) preferable from a computational standpoint. The particular case in which M depends 
linearly on p and polynomially on x, i.e., y = p~^ ip[x) where p{x) is a vector of monomials 
in the components of x, enables evaluating (21) and (22) analytically. In such a case, we 
obtain 


d{p, Si) 


\yi-p^^{xj)\ 

\\ph 


(24) 


d{p, Si,m) 


Vi - (p{xj)'^p 
p{\xi\)'^m 


(25) 


Once the norm d, and the value of the parameter 6 have been chosen, the uncertainty model 
fp{p), thus, the RPM in (18) are fully prescribed. 

The formulation below targets RPM with optimal characteristics. Two performance 
metrics used to quantify desirable traits of the prediction will be considered. One of them is 


1 . rPmax 

= / fp{P,^){M{xi,p) -yi)'^dp, 

i = l Jpmin 

1 

= V ^ - Vif) . (26) 

i=l 

where E[-\ and V[-] are the expectation and variance operators. This metric, which evaluates 
the dispersion of the random process Ry{x) about the observations, is computationally cheap 
to evaluate since it only requires evaluating means and variances. Another metric is 

N 

J{0,fp,z) = -Y[fM{xi,p{e)){yi), (27) 

i=l 


where fM(xi,p){yi) is the PDF of the predicted output evaluated oX y = yi. This density 
function depends on fp, which in turn, depends on 6. This performance metric is proportional 
to the probability of M reproducing the data sequence when fp[e){p) is the joint PDF of p 
and the observations in z are IID. Note that (27) is a likelihood function where 9, not p, is the 
independent variable. The sign in (27) enables maximizing the likelihood via minimization. 
The calculation of (27) requires constructing an accurate approximation of fy{y,Xi) for all 
Xj’s in z, a task that is more computationally expensive than the calculation of (26). The 
performance metric (26) is well suited to describe unimodal random processes. For instance, 
if z = {(xi,|/i), (xi,|/ 2 )}, the minimization of (26) will tend to yield a symmetric unimodal 
PDF for y{xi) centered about {yi + 1/2)/2. Conversely, the metric in (27) is better suited to 
capture multi-modal processes. For instance, if z = {{xi,yi), (xi,y 2 )}, the minimization of 
(27) will tend to yield a bimodal and symmetric PDF for y{xi) peaking at yi and p 2 - The 
modality of the predicted random variable y is exclusively dependent on the structure of M 
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and the spread of the observations. The framework above enables us posing an optimization 
program for calculating RPMs. 

Optimization Program 3. Consider the RPM defined by (18-20) for a norm d prescribed 
by any of the Equations (21-23) and any of the performance metrics in (26-27). Denote by 
0 the set of 6 values leading to a feasible P. The value of the parameter 9 required to fully 
prescribe this RPAfi is given by 

9 = argmin (j( 6 »,/p,z) : y{xi,P{9)) <yi< y{xi,P{9)), l<i<N] . (28) 

6» e 0 

Therefore, this formulation targets the support set of the PDF in (19) such that (i) 
the resulting random process Ry{x) maximizes the performance of the prediction, e.g., it 
minimizes the dispersion of probability about the observations, and (ii) all the observations 
in z can be reproduced by the computational model M{x,p) for at least one parameter 
point of P. Condition (ii), which is enforced via the set of inequality constraints in (28), is 
equivalent to P{9) fl S'* 7 ^ 0 for alH = 1, . . . , iV. Either of the parameterizations of P in (9) 
and (15) can be used in (28). In the context of (9) and (15), the constraint 9 E Q leads to 
a > 0 , and to pmin < Pmax respectively. 

2. RPM Example 

Next we calculate a RPM based on the same data sequence z^ and the same computational 
model M^{x^,p) used previously. Figure 6 shows the hgure of merit y(p) dehned by (20) and 
(23) for iV = 1 (top-left), N = 2 (top-right), N = 3 (bottom- left), and = 50 experimental 
observations (bottom-right). The top-left plot shows that the manifold Si, which is the set of 
parameter points where y(p) takes on the maximum value, is a nonlinear function of p. The 
rate of decay in the value of y(p), thus of fp{p), is driven by the rate at which the function 
{y\ — M^(p,x\)y‘ departs from its zero manifold. This rate can be adjusted via parameters 
fix and (82 ■ The top-right hgure is calculated by adding an additional term to the hgure of 
merit shown in the top-left hgure. This yields to a function y(p) where both manifolds Si 
and S2 are superimposed. The bottom-left hgure shows the hgure of merit corresponding 
to S'!, S2 and S3. Note that y(p) reaches its largest values where the zero-prediction-error 
manifolds intersect. The bottom-right plot shows the hgure of merit corresponding to all 
50 manifolds. Note that two-regions of high-manifold concentration emerge, with the one in 
the vicinity of p = [ 1 , 1 ]^ being dominant. 

Figure 7 illustrates the dependency of the hgure of merit y(p) on the tuning parameter 
(82 for N = 50. Values of (82 equal to one (top-left), 0.1 (top-right), 0.01 (bottom-right), 
and 0.001 (bottom-left) were used. Whereas larger values of (82 yield faster rates of decay 
and therefore, more accurate probability allocations (i.e., allocations that concentrate more 
probability closer to the manifolds), smaller values yield a smoother y(p), thus, density 
functions fp{p) that can be sampled more efficiently. Figures 6 and 7 are chosen to illustrate 
the (82 dependency of y(p) and are not illustrative of an optimal support set for p. 

A RPM dehned by Equations (18-20, 23, 26), with support set P given by (15) and (28) 
was calculated. The prediction corresponding to the resulting RPM, to be called RPM2, is 
shown Figure 8 . This formulation aims at concentrating the probability of the predicted re- 
sponse as close as possible to the observations. The comparison of Figures 5 and 8 illustrates 

^The parameters /3i and j32 in (20) used to prescribe fp(p) could be additional optimization variables. 
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Figure 6. Dependence of 7 (p) on the number of observations N. 



Figure 7. Dependence of 7 (p) on the tuning parameter /32- 
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that the mean response and the 5-percentile curves corresponding to RPM2 are a tighter rep- 
resentation of the data. Note that the mean prediction and the LS prediction are close even 
though they have a different mathematical structure (i.e., one corresponds to the evaluation 
of the computational model at a single parameter point whereas the other one is the average 
value of the model predictions corresponding to a large number of parameter points). The 
support set of p corresponding to RPM2, is bounded by Pmin = [0.8030, 0.8106]^ and 

Pmax = [1.0356, 1.2247]"'', for which J(9,fp,z^) = 3.8140. This region contains but it is not 
centered about points where fp{p) takes on its largest values. Conversely, the performance 
of RPMl is J{6, fp,z^) = 8.1332. This poor performance is expected since the probabilistic 
features of RPMl are not optimal. 


Validation Experiment 1 



Figure 8. RPM with minimal dispersion about observations (RPM2). 


V. Uncertainty Roll-up 

Uncertainty quantihcation aims at studying the effects of uncertainty in full-system tar- 
get applications. Unfortunately, we often have some tests at the subsystem/component level 
but not at the full system level. We want to use these tests to calibrate/ validate the compu- 
tational models of the subsystems in order to generate informative predictions of the target 
application. By roll-up we mean the process of calibrating the parameters of one or several 
computational models and using the resulting characterization of p to make predictions of a 
target application which also depends on p. In particular, we want to prescribe p according 
to the ability of the computational model M^{x^,p) to reproduce the data sequence for 
all j = 1, . . . ,n^. Note that the structure of the DGM and of the computational model 
corresponding to a validation experiment might be different from those for other validation 
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experiments. Once p is characterized, a prediction for the target application = M^{x^,p) 
will be made. Note that the response quantities y^, . . ., i/”” and y^; as well as the cor- 
responding states . . ., x"”' and x* may differ from each other (they might describe 

different physics and different quantities of interest). 

The IPM and RPM formulations above can be naturally extended to study the roll-up 
problem. The main extension considered in this section is the joint usage of all the data 
sequences with all computational models. The resulting characterization of p along with 
each computational model M* constitutes a predictor model. Therefore, we will obtain a 
family of predictor models linked by the common prescription of the uncertainty in p where, 
each family member describes one of the validation experiments. Extensions to the 
formulations above are presented next. 

Optimization Program 4. Consider the family of IPMs defined by (3-9) for M^(x^,p), 
M^{x‘^,p), . . M’^’'(x””,p) and assume that the values c and m > 0 are set in advance. The 
value of the parameter a required to fully prescribe this family is given by 

a = argmin | ^ P{9))] _ ,P{9)) < y(j < f (xh ,P(9)), (29) 

a>0 ny iJ-' - 

t J=1 

1 < j < ny, 1 < P < N^, where 6 = {c, m, a} } , 
and is a normalization constant. 

Therefore, (29) yields a P that minimizes the average spread of the Uy predicted responses, 
such that the envelopes of all such responses contain all observations. The normalization 
constant is used to make the spread of different response metrics comparable. One pos- 
sible choice for this constant is the sample variance of the observed responses corresponding 
to the jth validation experiment, = V[y^], where y{, . . . ,y(^j are the sample points. As 
expected, the spread of the prediction envelopes for all Uy models, as well as the size of P, 
grow with the value of Uy. The extension of Optimization Program 2, which follows the same 
structure and rationale of (29), can be easily inferred. Regarding RPMs, the extension of 
Optimization Program 3 is as follows. 

Optimization Program 5. Consider the family of RPMs defined by (18-20) for M^(x^,p), 
M^(x^,p), . . ., M"’'(x"”',p), and a norm d prescribed in any of the Equations (21-23). The 
value of the parameter 9 required to fully prescribe this family is given by 

0 = argmin {J(0,/p,z\...,z"’') : y^{xh,P{9)) < yh <y^{xh,P{9)), 
e e 0 

l<j<ny, 1<P<N^}, (30) 

where the cost function J is given by 



• • ’ 1 £)j ’ 

i=i 

(31) 

J{0,fp,z\ . 

TLy 

t=i 

(32) 
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when the performance metric J is given by (26) or (27) respectively. The parameters q > 
0 G and ?7 > 0 G M"” are constants set by the analyst, while is the normalization 
constant introduced earlier. 

Therefore, (30) yields a value of 6 for which the corresponding fp{p) minimizes the cost 
function J such that the envelopes of the predicted responses contain all corresponding 
observations. The cost J in Equation (31) is the weighted sum of the dispersion of all RPMs 
about the observations. Use g = 1 to give the same importance to the calibration of all 
RPMs. Conversely, the cost J in Equation (32) depends on the likelihood of all RPMs 
reproducing the corresponding data sequence. Use rjj = 1 to weigh each observation point 
equally, and pj = 1/N^ to weigh each data sequence equally. The spread of the envelopes 
and the tightness of the random predictions resulting from (30), regardless of the particular 
performance function used, can only degrade as the value of increases. As before, either 
of the parameterizations of P in (9) and (15) can be used in (30). In the context of (9) and 
(15), the constraint 0 G 0 leads to a > 0, and to P min < Pmax respectively. 

3. Roll-up Example 

We now consider Validation Experiment 2. The DGM associated to this experiment, from 
which we obtained the data sequence with = 50 observations, depends nonlinearly on 
the state and parameters, and it is subject to state-dependent randomness (see Appendix). 
Figure 9 shows the observations. The computational model M^(p, x“^), for which G M and 
p EMf, was built to describe this DGM. As compared to the DGM, the model is subject 
to epistemic-, aleatory- and model-form uncertainty. 

We will calculate a RPM based on and M^(p, x^) Erst. This RPM is defined by (18- 
20) with norm d given by (23) and performance metric given by (26). The corresponding 
LS parameter estimate, pls 2 = [1.3090, 1.0866]^, differs significantly from the LS parameter 
estimate for Validation Experiment 1, pi^si = [0.9420, 1.0689]^ (see Fixed-parameters IPM 
Example). This difference is caused by the uncertainty affecting the models and the noise 
affecting the DGMs. The resulting RPM, to be called RPM3, is shown in Figure 9. The 
performance of RPM3, J{9,fp,z^) = 77.6056, is much larger than that of RPM2. The 
support set of p corresponding to RPM3, P{p 2 ), is given by pmin = [1.1333, 0.3192]^ and 
Pmax = [1.5174, 1.8015]"'^. Note that the uncertainty affecting is signihcantly larger than 
that of M^. 

The simultaneous calibration of and using z^ and z^ is considered next. The 
mathematical structure of the DGMs associated to each validation experiment, thus of the 
data sequences, and of the corresponding computational models is significantly different. 
RPMs for both validation experiments were calculated using Equations (23, 26, 30) for p = 1. 
These RPMs, to be called RPM4 and RPM5, are shown in Figures 10 and 11. These hgures 
show the 5-percentile curves, the IPM envelopes, and the LS predictions corresponding to 
each validation experiment alone. The large differences between p^si and Pls 2 , and between 
the corresponding predictions illustrate that the uncertainty for the joint system is much 
larger than the uncertainty resulting from calibrating each computational model individually. 
The support set of both RPM4 and RPM5, P{0^), is bounded by pmin = [0.8049, 0.3413]^ 
and Pmax = [1.5173, 1.8018]"'^. Note that P{9‘i) ^ P{9\) U P{92), and the volume of P{9s) is 
much larger than the other two. 
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Output 


Validation Experiment 2 



Figure 9. RPM for Validation Experiment 2 (RPM3). 


Validation Experiment 1 



Figure 10. RPM for Validation Experiment 1 (RPM4). 
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Validation Experiment 2 



Figure 11. RPM for Validation Experiment 2 (RPM5). 


The comparison of Figures 8 and 10 indicates a considerable degradation in the quality 
of the prediction. The performance of RPM4, J{6,fp,z^) = 30.2080, is 7.92 times larger 
than RPM2’s; whereas the performance of RPM5, J{9,fp,z^) = 211.4523, is about 2.72 
time larger than RPM3’s. The degradation in the tightness of the prediction is caused by 
calibrating and jointly. 

Figure 12 provides insight into the process by which the probabilistic description of 
the uncertainty p is constructed. This hgure shows the joint PDFs corresponding to the 
RPMs resulting from calibrating alone (left-subplot), from calibrating alone (center- 
subplot), and from calibrating and together (right-subplot). The concentration of 
probability indicates how the zero-prediction error manifolds are distributed in the parameter 
space. 

Finally, we use the characterizations of p obtained above to make predictions on a target 
application. Figure 13 shows the RPM for the target application = M^{x^,p) when the 
uncertainty in p is characterized based on the hrst validation experiment (left subplot), 
and on the hrst and second validation experiments (right subplot). The LS predictions 
corresponding to the hrst (blue), and second (black) validation experiments, the 5-percentile 
curves (red-dotted), process support (red-dashed) and the expected output (red-solid) are 
shown. The degradation in the predictive response, in terms of the spread of the support 
and of the dispersion of probability, is apparent. Note that the predicted random process 
is skewed towards smaller output values. Further notice that whereas the spread of the 
output increases when both validation experiments are used for calibration, the variation of 
the mean output response has actually decreased. Notice that the prediction for the target 
application depends on the ability of to reproduce the sequence z^, on the ability of 
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Figure 13. RPM of target application for calibrations based on 1st (left) and both (right) 
validation experiments. Same line conventions used previously apply. 
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to reproduce the sequence z^, and on the structure of M*. 


VI. Conclusions 

This paper proposes a new paradigm for the calibration of computational models and a 
framework for uncertainty roll-up. The formulations proposed render a description of the 
uncertainty in the models parameters, thereby enabling the evaluation of the dispersion in 
model’s predictions. The process by which the predictor models are generated is based on 
modeling the discrepancy between model predictions and observations using deterministic 
and probabilistic means. In contrast to standard Bayesian approaches, the predictions re- 
sulting from IPMs and RPMs will not converge to a deterministic response function as more 
data is made available. Instead, IPMs converge to an interval valued function matching the 
spread of the experimental data, whereas RPMs converge to a random process matching 
their distribution. 


Appendix 

The structure of the DGMs used to generate synthetic data, and of the corresponding 
computational models are provided below. For validation experiment 1 we used 

y =x (pi -h exp(p 22 :)) + (1 - 0.5u)p2 cos(5pix) -h PiP 2 (x - tanh((pi P 2 )x) + 1)-F 
77(0.1 + 1.5|x|))(l-u), 

where x is uniformly distributed in X = [—1, 1], u = 0 for the DGM, u = 1 for M^, and rj is 
a standard normal. For the validation experiment 2 we used 

y =x (pip2exp(p2(3 1.2(1 - u))/4) - 2p2Cos(pi(4 1.2(1 - u))x)) + P 2 + 

PiP 2 (x + 1.2(1 - u) tanh(pi -h 3p2/4)a:), 

where x is uniformly distributed in X = [—2,1]. The DGM was created using u = 0 
and making p 2 a triangular random variable with support [0.75, 1.25] and mean value of 1, 
whereas was created using u = 1 for a scalar p 2 - Note that knowledge of the structure of 
the DGMs and of the models is not required, and for practical purposes they can be thought 
of as black boxes. 
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